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Abstract 



A Monte Carlo generator of the final state of hadrons emitted from an ultrarela- 
tivistic nuclear collision is introduced. An important feature of the generator is a 
possible fragmentation of the fireball and emission of the hadrons from fragments. 
Phase space distribution of the fragments is based on the blast wave model extended 
to azimuthally non-symmetric fireballs. Parameters of the model can be tuned and 
this allows to generate final states from various kinds of fireballs. A facultative out- 
put in the OSCAR1999A format allows for a comprehensive analysis of phase-space 
OO , distributions and/or use as an input for an afterburner. 
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Program summary URL: http://www.fpv.umb.sk/~tomasik/dragonl 
Program obtainable from: http://www.fpv.umb.sk/~tomasik/dragon 
RAM required to execute with typical data: 100 Mbytes 
Number of processors used: 1 

Computer (s) for which the program has been designed: PC Pentium 4, though 
no particular tuning for this machine was performed. 

Operating system(s) for which the program has been designed: Linux; the pro- 
gram has been successfully run on Gentoo Linux 2.6, RedHat Linux 9, Debian 
Linux 4.0, all with g++ compiler. It also ran successfully on MS Windows un- 
der Microsoft Visual C++ 2008 Express Edition as well as under cygwin/g++. 
Programming language: C++ 
Size of the package: 32 818 bytes 
Distribution format: tarred and gzipped archive 
Number of lines in distributed program, including test data etc.: 6368 
Number of bytes in distributed program, including test data etc.: 153 939 
Nature of physical problem: Deconfmed matter produced in ultrarelativistic 
nuclear collisions expands and cools down and eventually returns into the con- 
fined phase. If the expansion is fast, the fireball could fragment either due to 
spinodal decomposition or due to suddenly arising bulk viscous force. Particle 
abundances are reasonably well described with just a few parameters within 
the statistical approach. Momentum spectra integrated over many events can 
be interpreted as produced from an expanding and locally thermalised fireball. 
The present Monte Carlo model unifies these approaches: fireball decays into 
fragments of some characteristic size. The fragments recede from each other 
as given by the pre-existing expansion of the fireball. They subsequently emit 
stable and unstable hadrons with momenta generated according to thermal 
distribution. Resonances then decay and their daughters acquire momenta as 
dictated by decay kinematics. 

Method of solving the problem: The Monte Carlo generator repeats a loop in 
which it generates individual events. First, sizes of fragments are generated. 
Then the fragments are placed within the decaying fireball and their velocities 
are determined from the one-to-one correspondence between the position and 
the expansion velocity in the blast wave model. Since hadrons may be emitted 
from fragments as well as from the remaining bulk fireball, first those from 
the bulk are generated according to the blast wave model. Then, hadron pro- 
duction from the fragments is treated. Each hadron is generated in the rest 
frame of the fragment and then boosted to the global frame. Finally, after all 
directly produced hadrons are generated, resonance decay channels are chosen 
and the momenta and positions of final state hadrons are determined. 
Typical running time: Generation of 10 4 events can take anything between 2 
hours to a couple of days. This depends mainly on the size and density of frag- 
ments. Simulations with small fragments may be very slow. At the beginning 
of a run there is a period of up to 1 hour in which the program calculates 
thermal weights due to statistical model. This period is long if many species 
are included in the simulation. 
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1 Introduction 



In ultrarelativistic collisions of heavy atomic nuclei matter is probed at high 
temperature and density. The fireball thus created exists only for very short 
period of time; it quickly expands and hadrons decouple from it. If the collision 
energy is high enough, a phase of deconfined quarks and gluons is reached in 
the early phase of the collision. Due to pre-existing longitudinal movement of 
the incident nucleons and the inner pressure of the matter, it swiftly expands 
in longitudinal as well as transverse direction. At certain moment the energy 
density becomes too low to justify deconfinement and the matter hadronizes. 
From lattice QCD calculations we know that the change from deconfined to 
confined matter is rapid though smooth crossover in the region of the phase 
diagram at low baryochemical potential [I] and it becomes first order phase 
transition from certain value of \ib upwards [2"f3"f4"] . 

It is important for our discussion that the passage through this transition is 
very fast. In such a case, equilibrium description may be inapplicable and the 
phase transition would not proceed in the same way as probed on lattice. In- 
stead, considerable supercooling can occur. In the region of the phase diagram 
where first order phase transition is expected even spinodal may be reached by 
the system if the expansion rate is bigger than the nucleation rate of bubbles 
of the new phase [5|6] . Then, the fireball disintegrates into droplets similarly 
to spinodal fragmentation known in liquid/gas nuclear phase transition [7j. 

Seemingly, such a mechanism should not work in the region of the phase dia- 
gram with smooth crossover. Nevertheless, it has been argued that an abrupt 
rise of bulk viscosity at T c can suddenly make the fireball very stiff and if strong 
expansion is present in such a moment it can drive the system into fragmenta- 
tion [S]. Hydrodynamic expansion in nuclear collisions is possibly unstable [9], 
thus fragmentation at hadronisation phase transition seems likely scenario. 

On the other hand, hydrodynamically inspired parametrisations [T0lllllll2f 131114 
combined with thermal models [T5fl6|ll7irT8] provide satisfying description of 
single-particle spectra, particle abundances, and some also femtoscopy mea- 
surements. It is necessary to note that these are observables which are ex- 
tracted from data summed over a large number of events. Event-by-event fluc- 
tuations of mean p t [19,20,21] and angular correlations [22] indicate possible 
presence of clusters in momentum distributions. Such clusters could be due 
to fragmentation at the phase transition. Note that momentum clusters are 
buried under many entries to the histograms if summed over many events. 

It is the purpose of the present Monte Carlo droplet generator to produce 
artificial data sets which resemble those coming from real nuclear collisions 
provided fragmentation occurs at hadronisation and hadrons are emitted from 
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fragments without any further scattering. Its name, DRAGON, stands for 
DRoplet and hAdron Generator for Nuclear collisions. In a way, the model 
is similar to THERMINATOR [23J, with the crucial difference that emission 
from fragments is included. Note that the code can write out the final state in 
OSCAR1999A format [21] and thus a possible further evolution of the hadronic 
cloud can be simulated with the help of a cascade generator. 

In the next Section the model of particle emission from a fragmented fireball 
is reviewed. Section [3] explains the architecture of the program. Sections H] and 
[5] explain how to install and run the generator. Some representative results are 
presented in Section [6] and the paper in concluded in Section [7J Details about 
generation of momenta from Boltzmann distribution are summarised in the 
Appendix. 



2 The model of particle emission from fragmented fireball 



In rapid phase changes the fireball can fragment and hadrons are emitted from 
the produced fragments. There may also be a portion of the produced hadrons 
which is emitted directly from the bulk fireball. 



2. 1 Hadrons produced from bulk 



Directly produced hadrons are described by the blast-wave model. The Wigner 
distribution of their emission points and momenta is given as [T3"f2"5] 

S(x,p) d A x = m t cosh(y-r)) exp (-^^j 

x 6(1 — f(r, (j))) H(rj) 5(r — T )dr r drj r dr d(j) ■ (1) 

The model is formulated in terms of relativistic and polar coordinates r, 0, 77, 
and r 



x° = Tcosh?7 (2a) 

x 1 = rcos(p (2b) 

x 2 = rsin0 (2c) 

x 3 = t sinh r\ , (2d) 

while for momentum we use rapidity y, transverse mass (momentum) m t (pt), 
and azimuthal angle ip 
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p° = m t coshy (3a) 

p 1 =p t cos ip (3b) 

p 2 = p t sin ip (3c) 

p 3 = m t sinh y . (3d) 

Emission points are distributed uniformly in radial direction for 



where 

J? 

R x = aR, R y = — , (5) 

with i? being the mean transverse radius of the ellipsoidal fireball and a its 
spatial deformation parameter. 

The distribution H(rj) specifies the profile of the fireball in space-time rapidity. 
Wigner density in eq.([T]) is written in Boltzmann approximation and the factor 
p^u 11 gives the energy of the produced hadron in the rest frame of the moving 
fluid. The fluid velocity field is parametrised as 

= (cosh?? cosh^j, cos^b sinh^j, sin06 sinhr/t, sinhr/ coshr^) , (6) 
tan 0^ = a 4 tan , (7) 
77 t = fpov / 2(l + p 2 cos(20 fe )) . (8) 

Finally, the factor (2s + 1) in eq. (OQ) stands for the degeneration due to spin. 

Formulated in this way, the fireball can have elliptic transverse shape con- 
trolled by the parameter a and elliptic transverse flow profile parameterized 
by p2- 

2.2 Hadrons emitted by fragments 



The fireball decays into fragments of spherical shape (in their rest frame). They 
are placed according to the distribution (pQ) with T — > 0, i.e. their velocity is 
identical the local fluid velocity at the place at which they were produced. 

Fragments may come in one given volume b, if they are produced by mecha- 
nism which leads to one length scale which dominates the fragmentation — like 
spinodal fragmentation. Another possibility, implemented in the model is that 
the volumes are distributed according to gamma-distribution 

MV) = Jtt, (t)"' 1 exp (~) , (9) 



bT(k) \ b J "V b 
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with k = 2 [26pn . 



Fragments decay into hadrons exponentially in time, so the distribution of the 
emission time of a hadron in the rest frame of a fragment is 

^) = ^exp(-^) (10) 

where is the radius of the droplet (fragment). Hadrons are produced from 
the whole volume of the fragment with uniform probability. Their momentum 
is chosen according to Boltzmann distribution with the temperature in the 
rest frame of the fragment. 



2.3 Resonance decays 



Resonances may be produced from the fireball. Their lifetime is random ac- 
cording to exponential decay law exp(— Ft) (in the rest frame of the reso- 
nance) . 

If a resonance with mass M decays via two-body decay to daughters with 
masses mi and m2, their energies will be 



_ M 2 — mk + m? , , 

Ei = 2 — — t 11 

„ M 2 — m? + m 2 , , , 

E 2 = \— — 2 - 12 

and they will be receding back-to-back (in the rest-frame of the resonance) 
with momenta 

J{M 2 - (mi + m 2 ) 2 ) (M 2 - (mi - m 2 ) 2 ) 

|Pi| = N = ^ ^ • (13) 



In case of three-body decays, in the rest-frame of the resonance, all the mo- 
menta of the daughter particles lay within a plane. Under the assumption 
of transition amplitude independent of momenta, the energy distributions of 
daughter particles are uniform. Thus there is some freedom in choosing the 
energies and momenta of daughter particles in such a way that the energy and 
momentum are conserved 



Note the different use of the parameter b here and in [26]: while there it has the 
dimension of inverse volume, here it has the dimension of volume. This brings it 
here on equal footing with the case when all fragments have the same volume. 
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E x + E 2 + E 3 = M (14a) 
\Pi\ 2 + \P2 1 2 + 2 |pi| \p 2 \ cos6»i 2 = p 3 (14b) 

where Ei = \j\pi\ 2 + raf and #12 is the angle between the momenta pi and p 2 . 

Daughter particles from a resonance decay may also be unstable. In that case, 
they will decay according to the same procedure. 



2.4 Chemical composition 



Relative abundances of the individual species follow the prescription of chemi- 
cal equilibrium with temperature T c h and chemical potentials for baryon num- 
ber and strangeness and ^5. The density of species i is then given as 



rii{T ch , ij, b , us) = gi 



d 3 p 



cxp 



p 2 + mf - (fi B Bi + n s Si 



T, 



ch 



9i T 3 j ( m i _V± 

2^ ch vT cft ' T ch 



(15) 



where the upper (lower) sign is for bosons (fermions), is the degeneracy 
factor, and 



J-ch J-ch 



exp 



A 



x 2 + 



T 2 



T, 



ch 



-1 



IM = ^BBi + fXsSi ■ 

From this, the probability that a random particle belongs to species i is 

frr \ ni(Td» hb, ps) 
Wi{T ch , n B , ii s ) = — — , 

}2i rii(T ch , [i B , Us) 
where the sum in the denominator runs through all the species. 



(16) 
(17) 

(18) 



3 Programming structure and solution 



The structure of the program is outlined in Figured! It is written in C++. 
The code is split into three files: dropem . cpp includes the main function and 
functions for reading and initialising resonance decays and for their decay. The 
file dgener . cpp defines classes needed for the working of the MC generator, 
and specrel . cpp defines basic vectors and tensors and operations with them. 
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Fig. 1. The structure of the program. 



3.1 Introductory phase 



Parameters of the model and steering constants for compiling and running 
are specified in the file params.hpp. This also includes the list of all species 
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called pproperties. This is an array, whose entries are records for individual 
species. In the current version we standardly include all baryons with masses 
up to 2.0 GeV/c 2 and mesons up to 1.5 GeV/c 2 . 

Resonance decays are listed in a separate input file resonances . input. A 
file with this name must be present in the directory in which the program is 
running. It lists the modes in which resonances can decay. If at least one decay 
mode for given species is listed in this file, that species is treated as an unstable 
resonance and, vice versa, if species is not listed in the resonance file, it is 
stable. In order to store data of all decay modes of a resonance and simulate the 
decay of resonances, a class DecayPattern is defined in file dgener . cpp. After 
all decay prescriptions are read in, a link to the corresponding DecayPattern is 
added to each record of unstable species in pproperties. Vice versa, daughter 
particle species defined in the decay modes are linked to their properties in 
pproperties. 

Subsequently, weights for generation of all species are calculated according to 
eq. (TT8I) and stored for each sort of hadrons in the list pproperties. 

At the end of the initial phase, average energy of hadron at the specified 
chemical composition and kinetic temperature is determined. This allows to 
translate the expected multiplicity (specified as input parameter) into the 
energy which is contained in the fireball or the fragments. With the help of 
the energy density this is translated into the expected volume of the fragments. 
This later controls the number of generated fragments. 

3.2 Loop over events 

The number of events to be generated is given as a macro NOEvents in pa- 
rameter file params.hpp. The loop is repeated this number of times. 

DRAGON can generate fireballs with ellipsoidal cross-section, so first the di- 
rection of the event plane is chosen randomly. Thus each event has a different 
reaction plane as it is the case in real data. 

In the next step, sizes of fragments are either generated according to gamma- 
distribution, eq. or they are set all to the same value. Their number is 
chosen so that the energy which they contain corresponds to the multiplicity 
to be produced from fragments. 

The sizes, if they are not all equal, are then sorted with the help of a quick- 
sort algorithm. This is done in order to facilitate placing of the fragments in 
the reaction volume. Their positions are generated randomly as specified in 
Section 12.21 Rapidity of the fragments is chosen randomly according to either 
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Gaussian distribution with specified mean and width, or uniform distribution 
between specified minimum and maximum. Transverse position is chosen uni- 
formly within the ellipsis with radii R x and R y . The direction of R x is parallel 
to the event plane. Velocity of the fragment is then specified by eq. For 
each fragment, after its position is generated, it is checked if it does not overlap 
with any of the previously allocated fragments. If so, the position is generated 
anew. The algorithm starts with large fragments and continues to the smaller 
ones; otherwise there could appear problems with placing the big fragments 
at the end. 

Once all fragments are placed within the fireball, hadrons emitted from the 
bulk (not from fragments) are generated. This begins with determining the 
type of hadron according to the probabilities calculated in the initial stage. 
After this, the thermal component of hadron momentum is generated accord- 
ing to Boltzmann distribution with temperature T^. Rejection method is used; 
more details can be found in the Appendix. The direction of the momentum is 
random with isotropic distribution. The position from which hadron is emitted 
is chosen according to the emission function in eq. ([!]) with the same rapidity 
distribution as was used for fragments. The position corresponds to some value 
of the collective expansion velocity via eq. (jH])- The generated momentum is 
then boosted by this velocity. It is also checked that the generated particle is 
not produced within some of the fragments. If so, its position is rejected and 
generated again. The whole procedure is repeated until the expected number 
of hadrons from the bulk is generated. 

Hadrons emitted from fragments are generated in a similar way as those from 
the bulk. In the rest frame of the fragment, energy and momentum are gen- 
erated according to Boltzmann distribution and the position according to ho- 
mogeneous distribution. The times of emission are distributed exponentially, 
as seen from eq. ffTDT) . Finally, hadron position and momentum are boosted 
by the velocity of the fragment and the position is also shifted by the initial 
position of the fragment. After the generation of each hadron it is checked 
whether the energy contained in all hadrons from given fragment so far does 
not exceed the total energy of the fragment. If it does, no more hadrons are 
generated from that fragment. Thus energy and momentum are not strictly 
conserved in hadron emission from fragments. This is not bothering as long as 
the energy and momentum of the fragment are unknown in the experiment. 

The next step is to decay resonances. This is performed by void function 
DecayResonances. All hadrons are stored in a First-In-First-Out stack of par- 
ticle records. The record for each particle includes a link to its decay pre- 
scription; if the particle is stable then no link is present. The stack also keeps 
the number of particles. The algorithm counts stable particles. It takes one 
particle from the stack. If it is stable, the particle is put back into the stack 
and the counter of stable particles is increased by one. If it is unstable, then 
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params.hpp 



dropem 



Fig. 2. Dependences of the source files. 

decay is initiated. The counter of stable particles is reset to and daughter 
particles are put into the stack. This is repeated until the counter of stable 
particles reaches the number of particles in the stack. 

To simulate the decay of a resonance the algorithm first chooses the decay 
channel according to branching ratios. Momenta of daughter particles in the 
rest frame of the fragment are determined as explained is Section 12.31 These 
momenta are first rotated so that their orientation is distributed isotropically 
and then boosted by the resonance velocity. Then, a time in which the reso- 
nance decays is generated according to the exponential decay law. The position 
in which the resonance finds itself at that time is the initial position of the 
daughter particles. 

When all resonances are decayed, generated hadrons are written into the out- 
put file. 



4 Installing DRAGON 

The droplet generator DRAGON is distributed as a package with three C++ 
header files (specrel .hpp, dgener. hpp, params.hpp), three C++ source files 
(specrel . cpp, dgener . cpp, dropem. cpp), an input file with resonance decays 
(resonances . input) and Makefile for easy compiling on Linux. 

The dependences of individual source files in the distribution package are 
illustrated in Figure [2J Their content is as follows: 

specrel contains prototypes (.hpp) and definitions (.cpp) of objects for 
four-vectors and three-vectors together with operations on them, like addi- 
tion and multiplication with Minkowski metric, where applicable. Also ten- 
sors and their operations are introduced. There are functions for boosting 
four-vectors to other reference frames. Also an object Particle is defined 
which stores properties, momentum and emission position of a particle. 
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dgener ( . hpp and . cpp file) introduces specific technical tools which are 

used in the random generator, 
dropem. cpp contains the main routine with the code according to algorithm 

outlined in Section [3J 
params .hpp is the header file with settings and values for parameters which 

are used in the generation. This file can be edited in order to compile and 

run the generator with a different set of parameters, 
resonances . input is input data file which should be present (under this 

name) in the directory in which the generator is running. 

The generator should compile and run on any system with standard C++ 
compiler. It has been tested with g++ on Gentoo Linux 2.6, RedHat Linux 9, 
Debian Linux 4.0. It also ran successfully on MS Windows under Microsoft 
Visual C++ 2008 Express Edition as well as under cygwin/g++. 



5 Running DRAGON 

On Linux, Makefile makes sure that executing the make command will com- 
pile, link and prepare the executable file dragon.exe correctly. If make cannot 
be used, then one should follow the commands as they appear in the Makefile. 
Visual C++ by Microsoft will compile and link the code if all . cpp and .hpp 
files are included. 

The executable dragon.exe is run with one or two parameters like e.g. 

> dragon.exe outfile.out dropinfo.out 

The first parameter is the name of the output file, where generated particles 
will be written out by the program, in an output format which can be chosen. 
The second parameter, which is optional, stores information about the droplets 
which have been generated in each event. If the second parameter is not given, 
this information is not written out. If parameters are omitted, then the default 
output file DRAG0N_events . out is chosen. A file resonances . input must be 
present in the directory where the program is being run. 

5. 1 Input 

Settings for the program and parameters are chosen before the compilation in 
the file params. hpp. They are also explained in the comments within the file. 
Here, they are explained in more detail: 

NOEvents is the number of events which shall be generated 
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RANDOMDROPLETS is a logical constant; which is set either to 1, if the volumes 
of fragments are to be set randomly according to eq. (Q, or to if all 
fragment volumes should be equal 

OVERLAPFORBIDDEN is by default set to 1; it can be set to if the user wants 
to allow an overlap of fragments. The latter choice is unrealistic, but may 
be useful in some special studies, since forbidding fragment overlap leads to 
some anticorrelation of momenta when hadrons from different droplets are 
usually unlikely to have similar momenta. 

GAUSSIANJIAPIDITY is a logical constant which is set to 1 if the rapidity 
profile of the fireball should be Gaussian and to for uniform rapidity profile 

ACCEPTANCECUT declares which hadrons should be written out to the output 
file. This can be used to formulate simple conditions for detector acceptance. 
The macro which is defined here is a logical statement which is evaluated 
for each particle before it is written out. Following kinematic variables can 
be used: rapidity (defined as variable yrap), transverse momentum (pT), or 
azimuthal angle (Phi). The condition can be, e.g. 

#define ACCEPTANCECUT ( (yrap>-l . )&&(yrap<=l . ) ) 
i.e. hadrons with rapidity between -1 and 1 are recorded. 

WRITEOUT is an optional string which is written in the header part of the 
output file if it is in OSCAR1999A format [21]. By default, it is left empty. If 
a newline command '\n' appears, it should be followed by which starts 
a comment line in OSCAR1999A format. 

short FORMAT is an identifier of the format of the output file. Currently, four 
possible output formats are predefined. They are explained in the comments. 
Note that in the OSCAR1999A format two numbers in each line are added 
to the standard output. The last number is 1 if the hadron described in that 
line comes from resonance decay and otherwise. The next to last number 
identifies the fragment from which the hadron (or its parent resonance) 
stems; and is -1 for hadrons not originating from any fragment. 

BIGMASS is set to 10. by default. This is the parameter B described in the 
Appendix below equation ( 1A.3I) . Its setting requires some fine-tuning and 
we do not recommend to change it. 

Anr is set to 20. by default. This is the parameter A described in the Ap- 
pendix in eq. (1A.6|) . It is not recommended to change it. 

double f otemp is the kinetic freeze-out temperature in units of GeV. 

double Tch is the chemical freeze-out temperature in units of GeV. 

double mub is the baryochemical potential in GeV. 

double mus is the chemical potential for strangeness in GeV. 

double huen is the energy density within the fragments, in units of GeV.fm -1 

double minrap is the minimum rapidity of fragments or directly produced 
hadrons which will be generated. 

double maxrap is the maximum rapidity of fragments or directly produced 
hadron which will be generated, so fragments and direct hadrons are gen- 
erated with rapidities between minrap and maxrap. 
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double dNdy _total is the dN/dy of all hadrons which roughly should be 
generated. This variable is useful in case of uniform rapidity distribution 
and is only used for calculation of N.total (see next). If it is not used 
there, it can be commented out. 

double N_total is the total expected multiplicity in the interval between 
rapidities minrap and maxrap. For a uniform rapidity distribution it can be 
conveniently calculated as 

double N_total = dNdy_total * (maxrap-minrap) ; 
It can also be set equal to a given number. 

double DropletPart is a parameter between and 1 which determines the 
fraction of all hadrons that stem from the decay of fragments. It is 1 if all 
hadrons are produced from fragments and if they all are emitted directly 
from the bulk fireball. 

double rapcenter is only relevant for Gaussian rapidity distribution and 
gives the rapidity of the centre of the rapidity profile. It is irrelevant if 
uniform rapidity distribution is chosen. 

double rapwidth is the width of Gaussian rapidity distribution and is ir- 
relevant if uniform rapidity distribution is chosen. 

double rb is the radius, in fermi, of the fireball which decays into fragments 
and/or hadrons. It appears as R in eq. (j5J). 

double a_space is the spatial anisotropy parameter a as it appears in 
eq. ©. 

double tau is the parameter r from eq. ([!]), in units of fm/c. 
double etaf is the parameter po from eq. (jHJ). 
double rho2 is the parameter p2 from eq. (jHJ). 

double b is the parameter b of the gamma-distribution of fragment sizes in 
eq. (|2D, if RANDOMDROPLETS is set to 1. In case that RANDOMDROPLETS is set 
to 0, this gives the volume of fragments. 

long int seed is the initial seed for pseudo-random generator. If this is set 
to 0, the generator is initiated with machine time. 

int NOSpec is the number of species which can be generated in the simula- 
tion. 

PChem pproperties [] is a vector of structures PChem which store the 
records of properties of individual species. The vector has as many entries 
as specified by NOSpec. For one species, these properties must be given (in 
this ordering), 

(1) Monte Carlo ID number of species according to Particle Data Group 
[27] ; integer 

(2) mass in GeV/c 2 ; double 

(3) baryon number; integer 

(4) strangeness; integer 

(5) 1 if the species is boson or if it is fermion 

(6) spin degeneracy; integer 

(7) put 1. here, double; (this will be calculated later by the program) 

(8) put another 1. here, double; (this will be calculated later by the pro- 
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Fig. 3. An excerpt from the file resonances . input. 



gram) 

(9) put -1 in the last position, integer; this will be also determined by 
the program and will link the species to its decay prescriptions; it will 
remain -1 for stable particles. 

Data entries for one species are divided by commas. Record for one species 

should be input within curly brackets. 

Decays of resonances are listed in the file resonances . input. The structure of 
the records in that file is as follows (an example of the file resonances . input 
is shown in Figure [3]) : 

• A record of all decay modes of one resonance starts with a line with three 
numbers: MC code (identifier) of the resonance, its mass in GeV, and its 
width in GeV. If -1 is put in the position of the mass, the code automatically 
reads in the mass from pproperties [] . 

• Record of a two-body decay contains five numbers. First, there is branching 
ratio for the decay channel, which must be multiplied by the Clebsch-Gordan 
coefficient, if applicable. This is followed by MC code of the first daughter 
particle and its mass, and the same for the second daughter particle. Again, 
if -1 is put instead of the mass, the mass it fed in automatically. 

• Three-body decays are recorded with seven numbers, where MC code of the 
third daughter particle and its mass are added to the structure of the record 
of two-body decays. 
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If the sum of branching ratios is not unity, the program will multiply them 
with a common factor so that their sum will become 1. Decays into photons 
(or leptons) can also be included with the appropriate MC codes (e.g. 22 
for photons). Any line that starts with or is empty is considered as a 
comment. Comments are very useful in order to enhance the readability of 
the records. 

5.2 Output 

Output data are directed into the file which is specified as a command line 
parameter. The possible formats for the output are explained as comments 
in params.hpp. Note that in case of the OSCAR1999A output format each 
line contains also two additional output numbers: On the 13th position the 
number of the fragment from which the hadron originates, or the number from 
which its parent resonance was emitted. For hadrons emitted from bulk this 
number is -1. On 14th position there is if this is directly produced hadron 
or 1 if the hadrons stems from a resonance decay. 

If droplet information is stored, the file will contain lines with the following 
structure: event number, number of droplet, rapidity of droplet, transverse 
velocity of droplet, azimuthal angle of its transverse velocity. This is followed 
by the multiplicities of individual sorts of particles which are emitted from the 
given droplet. Thus for example a record like 

# eid d_id rapidity v_t phi sorts: 

# -211 111 211 

#2 3 1.2724 0.3411 2.4513 12 10 16 

means that in event 2 we have droplet number 3, which moved with rapidity 
1.2724 and transverse velocity 0.3411c under azimuthal angle 2.4513 rad. From 
this droplet 12 7r~s (code -211), 10 tt°s (code 111), and 16 7r + s (code 211) are 
emitted. (Usually, one would make simulation with many more species; we 
only show three here for brevity.) 



6 Sample results 

6. 1 Size of droplets 

In order to get an impression of the size of droplets, we first give average 
numbers of final state hadrons coming from one droplet for various settings of 
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parameters. Resonances were included in this calculation. They decayed into 
stable hadrons and here the final numbers of stable hadrons are given. They 
were all calculated for energy density and chemical freeze-out conditions at 
RHIC T ch = 168 MeV, fi B = 266 MeV, fi s = 71 MeV. For small droplets with 
the volume of 2 fm 3 and kinetic temperature 160 MeV, there are on average 2.3 
pions, 0.2 nucleons, and 3.2 hadrons in total produced. For the volume 10 fm 3 
the corresponding numbers are 8.7, 0.9, and 12.1. For very large volume like 
100 fm 3 we have on average 80.9 pions, 7.8 nucleons, and 112.3 hadrons. The 
temperature dependence is rather weak, e.g. in the latter case dropping the 
temperature from 160 to 100 MeV increases the pion number to 94. In general, 
lower kinetic temperature leads to larger number of hadrons since less energy 
is used in the form of kinetic energy. 



6.2 Comparison with THERMINATOR 



DRAGON is similar to THERMINATOR [23J, with some differences: 

• DRAGON allows for particle emission from fragments; if fact this was the 
main motivation for conceiving it. 

• The radial profile of transverse expansion velocity grows linearly in DRAGON 
(cf. eq. (jSD), while it is kept at constant value if the blast wave model option 
is set in THERMINATOR. 

• THERMINATOR also simulates freeze-out in Cracow single freeze-out model 
[13] and in blast wave model with varying time dependences of radial coor- 
dinate of the freeze-out hypersurface [28J. 

• DRAGON allows to simulate azimuthally non-symmetric fireballs. 

• THERMINATOR uses SHARE [151129] to define chemical composition and 
resonance decays. The list of included resonances and decays is longer than 
that of DRAGON, which should not cause a problem, however, as higher 
lying states are suppressed by Boltzmannian factor. 

The results of both models are compared in Figure HI For this comparison, 
THERMINATOR was run with the set of parameters with which it is dis- 
tributed, just the freeze-out model was set to be blast wave and not the Cra- 
cow single freeze-out. Parameters of DRAGON were chosen correspondingly. 
Note that multiplicity results from the chosen parameters in THERMINATOR 
while it is set by hand in DRAGON. Thus the difference in absolute normal- 
isation of the spectra is irrelevant. Recall, that the transverse flow expansion 
rapidity 0.55 was constant irrespective of radial coordinate in THERMINA- 
TOR; in DRAGON we chose po = 0.55. This leads to the same mean transverse 
expansion velocity. 
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Fig. 4. Comparison of pt spectra generated by DRAGON (solid lines) and THER- 
MINATOR [23] (dotted lines). In THERMINATOR the blast- wave source was cho- 
sen with the parameters that come with the standard distribution: v r = 0.55, 
r = 9.74 fm/c, p max = 7.74 fm, T = 165.6 MeV, \i B = 28.5 MeV, fi s = 6.9 MeV. 
Parameters were chosen correspondingly in DRAGON with po = 0.55. Spectra from 
10 simulated events are shown. 

In Figure H] we observe this difference. At low p t , THERMINATOR spectra 
are suppressed against DRAGON. These are the hadrons with low transverse 
velocities. In DRAGON, they are produced mainly from regions of the fireball 
which move slowly outwards. In THERMINATOR, such regions are missing 
(since transverse velocity is everywhere the same) and this leads to the ef- 
fect on the spectra. Pions move relativistically already at rather low p t , so 
for them this suppression is at very low p t and is almost invisible. On the 
other hand, at very high transverse momentum their spectra become flatter 
in THERMINATOR. These are pions which move with very high transverse 
velocities — higher than the transverse expansion velocity of the fireball. In 
THERMINATOR, their production is enhanced by the fact that we have 
larger region that moves outwards with rather large transverse velocity. In 
other kinematic regions — high p t for kaons and protons and semi-low p t for 
pions — the spectra from DRAGON and THERMINATOR are parallel and 
consistent. 



6.3 Single-particle spectra 



Results obtained from a simulation with Gaussian rapidity profile, are illus- 
trated in Figures O O and [7J First of them shows the rapidity spectra. In 
the simulation, space-time rapidity distribution was centered at rj = with 
a width of 1.3 and kinetic freeze-out temperature Tk = 160 MeV. We show 
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Fig. 5. Comparison of rapidity spectra from a model with fragments (solid lines, 
b = 50 fm 3 ) and without them (dotted lines) for charged pions (top), kaons (middle), 
and protons and antiprotons (bottom). A model with Gaussian space-time profile 
was chosen with the width Arj = 1.3. Other parameters were: = 160 MeV, 
T ch = 168 MeV, \i B = 266 MeV, fi s = 71 MeV, p = 0.6, Spectra from 10 4 
simulated events are shown. 

results from simulations with and without fragments. Note, that in case of 
fragments their typical size is actually chosen rather large (b = 50 fm 3 ). No 
difference is observed between the spectra obtained from the two kinds of sim- 
ulation. This is expected: fragmentation of the fireball cannot be seen from 
observables integrated over large number of events. 
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Fig. 6. The pt spectra from a model with Gaussian space-time rapidity profile with 
the width Aw = 1.3. Results from the same simulation as in Fig. [5j = 160 MeV, 
T ch = 168 MeV, p B = 266 MeV, p s = 71 MeV, p = 0.6. Spectra from 10 4 simulated 
events are shown; they are integrated over rapidity. Solid lines show results from 
a simulation with fragments with b = 50 fm 3 and dotted lines show results from a 
simulation without fragments. 



The same conclusion can be drawn from the pt spectra. The spectra from non- 
fragmented fireball appear slightly natter. This is because in such a simulation 
particles can be produced from regions close to the radial edge of the fireball 
which move with highest transverse velocity. Fragments are produced in such 
a way that their volume is always completely included in the volume of the 
fireball and their velocity corresponds to the fireball expansion velocity at the 
centre of the fragment. Thus fragments (especially the big ones) can never 
obtain as high transverse velocity as the hadrons in non-fragmented fireball. 
In a case where both simulations are performed with the same po, like in Fig. [61 
simulation without fragmentation thus yields slightly flatter p t spectra. 

In Fig. [7| we show the contributions to pion spectra. In both cases, with and 
without the fragments, they are the same. The most important contribution 
at this chemical freeze-out temperature (T ch = 168 MeV, ps = 266 MeV, 
Ps = 71 MeV) comes from resonance decays. Direct pion production becomes 
equally important at around p t = 1.5 GeV/c. 
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Fig. 7. The contributions to pion spectra from the simulations of Fig. Thin solid 
line is the resonance contribution, dotted line shows direct pion production and 
thick solid line is the total resulting spectrum. Left: simulation without fragments 
(droplets), right: simulation with fragments (droplets) with b = 50 fm 3 . 

6.4 An illustration of fragments 



As mentioned, fragmentation cannot be observed in spectra integrated over a 
large number of events. Correlation and fluctuation observables go beyond the 
scope of this paper and shall be investigated in subsequent dedicated papers. 
Just to illustrate the production from fragments, momenta and positions of 
particles are illustrated in Figure [HJ Positions of emission points (right column) 
and momenta of hadrons (left column) in rapidity and azimuthal angle are 
shown. In this illustration only directly emitted pions have been considered 
for the sake of simplicity. In the top row we show momenta and positions of the 
emission for a fireball that expands boost invariantly in longitudinal direction. 
It is observed that both positions and momenta spread uniformly over the 
whole (77, 0) interval. In order to clearly illustrate the effect of fragments on the 
momentum distribution, in the bottom row results from an event with large 
fragment volume (b = 50 fm 3 ) and very low kinetic freeze-out temperature 
(Tfc = 10 MeV) are shown. The large volume ensures that there are many 
pions coming from one fragment and the low temperature puts a limit on their 
thermal velocity so that it stays close to that of the fragment. Thus clustering 
of the emission points is transferred to the momentum space. Note, however, 
that the values of parameters are unrealistic. In order to show results in a 
more realistic simulation, in the middle row a temperature = 170 MeV was 
chosen. The clustering in momentum space is still observed though now it is 
much less pronounced. Note that in realistic situation (i) the size of fragments 
will probably be smaller; (ii) the number of fragments will be bigger; (hi) also 
other hadron species shall be present and pions shall dominantly be produced 
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Fig. 8. Positions in space-time rapidity r] and azimuthal angle from which pions were 
emitted (right column) as well as rapidities and azimuthal angles of the momenta 
which they had (left column). No resonances and no other particles than pions were 
produced here. Top row: pions emitted from boost-invariantly expanding fireball 
which does not fragment. Middle row: pions emitted from two large fragments with 
temperature 170 MeV. Hadrons from the same fragment are indicated by the same 
symbol. Rapidities, transverse velocities, and azimuthal angles of the fragments are 
(-0.55, 0.51c, 3.53 rad) and (0.84, 0.47c, 3.30 rad) and they emit 139 and 123 pions. 
Bottom row: as middle row but for the sake of illustration the temperature was 
set to 10 MeV. Rapidities, transverse velocities and azimuthal angles of the two 
fragments are (-0.44, 0.45c, 4.49 rad) and (0.73, 0.34c, 0.92 rad) and they emit 435 
and 390 pions. 
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p ( [GeV/c] 

Fig. 9. The elliptic flow parameter v 2 of all hadrons calculated for parameters of the 
blast wave model that have been fitted to Au+Au collisions at y/s = 130 ^4GeV, cen- 
trality class for elliptic flow 11-45% [3D]: T k = 107 MeV, p = 0.85/V% P2 = 0.058, 
R = 11.1 fm, a = 0.939, r = 7.7 fm/c p3] and T ch = 174 MeV, ll b = 46 MeV, 
Us = 13.6 MeV |17j . Solid red lines with triangles show results from simulation 
with fragments (b = 10 fm 3 ) and dashed blue lines with squares show results from 
simulation where all hadrons have been produced from the bulk fireball. 



from resonance decays. Thus clustering in momentum space shall be even 
less visible and sophisticated methods might be necessary in order to identify 
fragmentation. 



6. 5 Elliptic flow 



An important observable is the elliptic flow. This model provides the possi- 
bility to simulate azimuthally non-symmetric collisions where elliptic flow is 
observed. Within the blast wave model, the existence of azimuthal anisotropy 
of single particle spectra can be achieved by setting the spatial anisotropy a 
and/or the expansion velocity anisotropy p 2 [25]. In order to illustrate the ellip- 
tic flow, we calculated v 2 for a set of parameters which reproduce mid-central 
Au+Au collisions at = 130 AGeV [H] (Figure EJ). 

The elliptic flow parameter t> 2 was determined from simulated data via two- 
particle correlations in the azimuthal angle 
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JdA cos(2A) N 2 (A, Pt ) 



1/2 



fdAN 2 {A,p t ) 



(19) 
(20) 



where N\(<p,pt) is the single-particle distribution in azimuthal angle and trans- 
verse momentum. In practice, v 2 is not evaluated for a sharp value of the 
transverse momentum, but momenta from an interval in p t are taken. In this 
analysis, the p t range was divided into seven intervals. 

In Fig. Owe compare the elliptic flow for all hadrons calculated in a standard 
blast wave model with one calculated in case of all hadrons emitted from 
fragments. In both cases, all parameters of the model were the same; only 
the percentage of particles coming from fragmets is in one and 100% in the 
other case. A slight increase of the elliptic flow in case of fragmented fireball is 
observed. This is due to enhanced correlations of hadrons from one fragment. 
A thorough investigation of the elliptic flow from fragmented fireball shall be 
deferred to a dedicated paper. 



7 Conclusions 



The new Monte Carlo generator of the final positions and momenta of hadrons 
represents the blast wave model and incorporates important additional fea- 
tures. Above all it is the possibility to simulate hadron emission from a frag- 
mented fireball. It is possible to vary the planned multiplicity and the per- 
centage of particles which are emitted from the fragments. Additionally, it 
also allows for simulation of azimuthally non-symmetric events. THERMINA- 
TOR, the closest related model on the market, exists only in the azimuthally 
symmetric version. 

Thus two kinds of studies can be envisaged for which DRAGON appears 
unique: (i) studies of correlations and fluctuations connected with fragmen- 
tation of the fireball; (ii) studies of spectra and correlations in non-central 
collisions. The main driving motivation for its development was the former, as 
fragmentation may be intimately connected with the dynamics of the phase 
transition. The Monte Carlo generator can be used for designing and testing 
new observables for probing fragmentation. 
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A Momentum generation from Boltzmann distribution 



In this appendix it is reviewed how the momenta of hadrons are generated 
according to Boltzmann distribution, in the rest-frame of the fragment or 
fluid element. 

First, the case = is treated separately. In such a case the momentum 
vanishes and the energy is put equal to the hadron mass. 

For non-zero temperatures, rejection method is used to generate the size of 
the momentum. The Boltzmann distribution is 

V B (p)ocp*e W (- Vma r +p2 ) . (A.l) 

If the mass of the hadron is small with respect to T, this distribution re- 
ceives large contribution from the region of high momenta. Momentum is first 
generated according to gamma distribution 



r 3 (p)ocp 2 exp(-|) , (A.2) 



which is implemented as explained in [21]. Once momentum p is generated, it 
is accepted with probability 

(p- y/rri 2 +p>\ 
V acc = exp I I . (A.3) 

The natural scale for the mass is the temperature. In the code, this procedure 
is used for masses smaller than B ■ T, where the parameter B is introduced as 
macro BIGMASS and is normally put equal to 10. Acceptation probabilities of 
some light particles are shown in Figure IA.1I 

For large masses, this procedure becomes ineffective because the acceptation 
probability is tiny if p <C m. In this case, the problem can be treated non- 
relativistically. All three components of the momentum are generated accord- 
ing to the non-relativistic distribution 

P m -M cx exp • (A.4) 

Thus total momentum will satisfy the distribution 

V nv (p) ccp 2 exp (^-^-^j . (A.5) 

This cannot be used in the rejection method in a strict way, because distribu- 
tion (1A.5I) drops for large p as e~ ap whereas the correct relativistic distribution 
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Fig. A.l. Top: The acceptance probability for light hadrons for some combinations 
of mass and temperature calculated from eq. (|A.3|) . Bottom: Boltzmann distribution 
of momenta for light hadrons. 

( 1A.1I) goes like e~ /3p . Thus there will always be a region of large p in which 
kV qt (p) becomes smaller than Vb{p) no matter the value of the multiplier k. 
The acceptance probability calculated as Vb(p)/ nV nr (p) then becomes bigger 
than 1. 

In the implementation, the constant k is chosen in such a way that the accep- 
tance probability becomes 1 at p = AT. Thus the acceptance probability then 
is 

This technically means that for p larger than AT all hadrons with momenta 
generated according to non-relativistic distribution are accepted. It also means 
that in the spectrum these are suppressed relative to the case where they would 
be generated according to the correct distribution, due to faster decrease of the 
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Fig. A. 2. Top: The acceptance probability for heavy hadrons for some combinations 
of mass and temperature calculated according to eq. (|A.6jl . Bottom: Boltzmann 
distribution of momenta for heavy hadrons. 



non-relativistic distribution. The parameter A must be chosen. If it is small, 
then as a result the algorithm will produce too few hadrons with momenta 
above AT. On the other hand, if A is too large, then the acceptance proba- 
bility will become very small at small momenta and the algorithm becomes 
ineffective or even not working. For usual freeze-out temperatures, between 
say 100 and 180 MeV, the default reasonable setting is A = 20. For this 
choice, in case T = 0.1 GeV and m = 3 GeV/c the fraction of hadrons which 
is suppressed due to this artifact is at the level 0.86%. Increasing tempera- 
ture and decreasing the mass makes this fraction yet smaller. On the other 
hand if this algorithm is used for masses larger than 10T, then minimum ac- 
ceptance probability is reached for the smallest mass 10T and p = and is 
about 4.8 • 10~ 4 . In Figure IA.2I the acceptance probabilities as functions of 
momentum for heavy (non-relativistic) hadrons are shown. 
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